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Population dynamics in random ecological networks are investigated by ana- 
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lyzing a simple deterministic equation. It is found that a sequence of abrupt 
changes of populations punctuating quiescent states characterize the long time 
behavior. An asymptotic analysis is developed by introducing a log-scaled 
time, and it is shown that such a dynamical process behaves as non-steady 
weak chaos in which population disturbances grow algebraically in time. Also, 



^ ■ some relevance of our study to taxonomic data of biological extinction is men- 

tioned. 
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Almost all possible equilibrium populations in random ecological networks become un- 
stable when the number of fully connected species is sufficiently large This is a general 
result obtained from random matrix theory . The loss of stability of equilibrium popula- 
tions, of course, dose not necessarily imply a catastrophe of the system; it is known that limit 
cycles, heteroclinic cycles, and chaos arise beyond the loss of the stability in a few species 
system [§-0]. Therefore, it seems natural to consider how populations evolve in ecological 
systems consisting of many randomly interacting species. 

This problem is related to the co-evolution of many species in an unchanging physical 
environment. Each species evolves in the effective environment formed by the others, and 
the evolution of one species alters the environment of many others. This leads to a kind 
of frustration in the ensemble of species, and if the frustration is never resolved by the 
evolution, it would continue indefinitely. Such a scenario was proposed by Van Valen || and 
is often referred as the Red Queen hypothesis [[7].|| . 

Our aim is to present significant features of population dynamics in a random ecological 
network by analyzing suitable mathematical models. The choice of models may be a con- 
troversial step, and various types of dynamical behavior will be observed depending on this 
choice. We postpone discussion about this step, and in this paper we focus on a specific 
model which is composed of simple deterministic equations for the populations of N species. 
The equations are so simple that we can employ an asymptotic analysis for a long time 
behavior of the populations. 

Our equations describe the time development of populations {xj} of N species. We 
assume that each species grows at a rate Aj which is expressed by a Lotka-Volterra type 
coupling with the other species: 

\i = cii + ^gijXj, (1) 
j 

where aj and gij correspond to an intrinsic growth rate and competing coefficients, respec- 
tively. We further assume that each population has a saturation level, and we normalize the 
populations so that Xi — 1 at the saturation levels. Xi then obeys the logistic equation: 
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dxi 
~dt 



\i{Xi Xj), 



(2) 



where < x, < 1, and x, = and Xj = 1 correspond to 'extinction' and 'saturation', 
respectively. Although we consider the case that a, = in (|l|) for simplicity, the analysis 
given below is easily extended so as to include this term. From our motivation, N is a given 
large number, and g^ is assumed to be a random variable obeying a Gaussian distribution 



renormalized to 1 by changing the time scale. 

Equation (Q) has 2 N equilibrium states, each of which is expressed by (xi, X2, • • • , xat) = 
(o"i, cr 2 , ■ • ■ crjv), where <jj G {0, 1}. In an equilibrium state (cr 1 , cr 2 , • ■ ■ , a at), species satisfying 
Oi = 1 exist at their saturation levels, while other species are all extinct. We call the special 
state with Xi = for all i 'perfect extinction'. The linear stability of equilibrium states 
is determined by the eigenvalues of the linearized equation around these states. It can be 
easily shown that the probability of choosing a stable state from all equilibrium states equals 
1/2^, and therefore almost all equilibrium states correspond to saddle points. Although the 
expectation of the number of stable equilibrium states is approximately 1, whether or not 
a given system can reach such a state depends on the global structure of the flow in phase 
space. 

We first give results of numerical simulations of ([]]) and (|2|) in Fig. |1| It is seen that a 
population abruptly increases or decreases from a quiescent state close to an equilibrium. 
We call the quiescent state and abrupt change 'quasi-equilibrium' and 'burst', respectively. 
Such dynamical behavior seems to continue indefinitely within our computationally available 
time. In fact, even if the system eventually settles into a stable equilibrium, a long transient 
time is needed to reach the neighborhood of this state. We can estimate the waiting time 
as 0[exp(2 JV )], because the system needs to experience 0(2^) bursts to meet the stable 
equilibrium state, and the time interval between two successive bursts grows exponentially 
on the average. (The latter statement will be demonstrated in the argument below.) We 
thus focus on the long time behavior composed of 'quasi-equilibrium' and 'burst' irrespective 



with zero mean and deviation J/ \/N except for g, 



u 



0. Note that the value of J can be 
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of the final state. 

Now, we develop a new asymptotic method. Our key idea is to introduce a new variable 
Ui defined as to satisfy 

Xi = f{yi) := ^ xp( ^) (3) 
Then, Equation (^|) is transformed into a simple form: 

< 4 » 

Integrating this equation, we obtain 

Vi(t) = yi(0)+ <\>t, (5) 

where < Aj > is the average of the growth rate over the time interval from to t and is 
written as 

<A,>=i f dt'X^f). (6) 
t Jo 

Since < Aj > does not converges to zero for t — > oo except in the case that the system 
approaches the perfect extinction, we obtain the asymptotic form of £j for t — > oo: 

^(*) = /(%(0)+ < Xi > t) -> 0(< A, >). (7) 

That is, after a sufficiently long time, the system stays at a quasi-equilibrium state, and 
whether Xi is close to or to 1 is determined by the sign of < Aj >. We say, in this 
sense, that the dynamics of the population is slaved to that of the average growth rate. We 
thus consider time development of < Aj >. The definition of < Aj > given by (|6|) leads 
immediately to the equation: 

j t = - t (~ < \ > +Y,9ijXj), (8) 

where we have used the expression of the growth rate ([!]). Then, by introducing a log-scaled 
time t = log(t), defining a new variable hi{r) =< Aj >, and using the asymptotic form ([7]), 
we obtain the equation describing the slow dynamics of the average growth rate: 



dhi 
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(9) 



Here, it is worth noting that the equations of motion for the average growth rates are 
autonomous when we use a log-scaled time r. We now consider the relationship between the 
orbits of hiij) and Xi(t). First, from the asymptotic form (|7|), x-i is close to or 1 depending 
on the sign of hi, and a burst of the i-th species is identified with a zero-crossing of hi. Second, 
consider periodic motion of hi expressed by hiij) = sin(27rr/T). Then, a log-scaled time at 
which hi crosses zero for the n-th time is given by r n = Tn/2. Correspondingly, zeros of 
< Xi > occur at exp(Tn/2). This means that the time interval between two successive bursts 
grows exponentially in n. We can see in general that a periodic orbit of hi(r) corresponds to 
a transient orbit of Xi(t) attracting to a heteroclinic cycle. Finally, the long time behavior 
of populations corresponding to chaotic motion of hi is described by a non-steady process 
of an irregular occurrence of bursts. The time interval between two successive bursts grows 
irregularly in n, but on the average grows exponentially. 

We also note that Equation (|9|) has a form similar to the random neural network model 
proposed by Sompolinsky et.al. which is obtained when we replace 9(hj) by tax\h(Khj). 
They studied the statistical properties of fluctuations by developing a new method called 
'dynamical mean field theory' and showed that their equation exhibits chaotic behavior when 
KJ > 1 in the limit N — > oo 0. Therefore, one may guess that ([|) has chaotic solutions 
for sufficiently large N. 

In order to investigate dynamical properties of solutions hi(r), we performed numerical 
simulations of (|S]). Noting that the equation is piece- wise linear and discontinuous at hi = 0, 
we can construct a solution starting from an initial condition {/ij(0)} in the following way. 
First, the equation @ can be integrated analytically until a time T\ at which the first burst 
occurs, and /ij(rx) is expressed by 



where Xi = J2j 9ij6(hj) takes a constant value during this period, and T\ is given by 



hi{n) = (hi(0) - Aj)exp(-ri) + A 



(10) 
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n = Min iog(i - ^®), (11) 

where the index i runs over the species satisfying hi(0)Xi < 0. By repeating a similar 
procedure, we obtain a sequence of r, with (i = 1, 2, • • •), and simultaneously the solution of 

(D- 

In Fig. 0, we plotted the species which bursts at r, with (z = 1,2, •••). The plotted 
pattern seems to be random along the time axis, but certain inhomogeneity in species can 
be seen. In fact, Fig. |3| shows that the time interval between two successive bursts obeys 
a Poisson distribution P(r), while as shown in Fig. £|, Q(r), the distribution of persistence 
time of positive hi, has a power law tail 

Q(r)~r- 2 (12) 

for t — > oo, which is due to the randomness of the long-time average of hi. Therefore, in 
the random neural network model proposed by Sompolinsky et.al. ||, similar behavior is 
never observed |fLl}| . From the forms of these distributions, we can expect that the dynamics 
of hi(r) are steady. As a result, on the real time scale, the population dynamics are non- 
steady; the time interval between two successive bursts grows in n irregularly, but on the 
average exponentially. Further we numerically calculated the maximum Lyapunov number 
/j, for many set gij chosen randomly, and found this number to be positivein each case. Thus, 
a disturbance 8hi{r) grows as 5hi(r) ~ exp(/ir), that is, time series of h(r) are orbitally 
unstable. This result leads to power law divergence of a disturbance of populations 6xi(t) 
of the form 

<M*)~*". ( 13 ) 

and thus zero Lyapunov exponents. Our population dynamics therefore are not chaotic, but 
'weakly chaotic' [K|. Discussion of the distribution function of \i for the ensemble of {gij}, 
the Lyapunov spectrum, their dependence on N, and a proof of (0) will be presented in a 



separate paper | XT 



We have shown the significant features of the long time behavior of populations described 
by ([!]) and (^]). These are expected to be common to a large class of models of random 
ecological networks. Recently, population dynamics in random ecological networks have 
been studied by analyzing another type of equation fl2| . |T3|| . This equation is obtained by 



normalizing the sum of the populations to unity, and the resultant form is the same as the 
replicator model proposed in the context of molecular evolution We have found that a 



non-steady process with a log-scaled time appears generally in a random network |14|, and 



one of the authors has investigated a transition to 'chaos with a log-scaled time' in a four 



dimensional system [15 



Here, we present a remark: a small perturbation of the dynamical system, e.g., the addi- 
tion of eF(xi) to the right-handed side of (|J), yields topologically unequivalent orbits Xi(t). 
In other words, the above solutions are structurally unstable. This structural instability 
is analogous to that of homoclinic orbits |16| |. Since structurally unstable behavior cannot 



be observed without certain constraints, one may object that the behavior we study is not 
generic. We believe however that an understanding of the behavior of the unperturbed 
system will provide the first step toward understanding the rich variety of dynamical phe- 
nomena of perturbative systems much as center manifold theory has succeeded in describing 
complex behavior near bifurcation points . 



We finally discuss the relevance of our study to the taxonomic data of biological extinc- 
tion. According to Van Valen |8j , the lifetime of species obeys a Poisson distribution with an 
extinction rate Q which is about 10 -7 [year -1 ]. (Such a simple picture has been called 'con- 
tinuous extinction'. The alternative to this is the 'episodic extinction' picture ||17|| . ) Since 
extinction (xi = 0) never occurs in our model, we need to introduce another assumption in 
order to see the correspondence with extinction in taxonomic data. One plausible idea is to 
regard an abrupt decrease of a population from near the saturation level as extinction. This 
may be justified, because the taxonomic data are made from fossil records. Then, noting 
that the time spent near the saturation level measured with a log-scaled time corresponds 
to the persistence time of positive hi, we find that the lifetime of a species obeys the dis- 



tribution Q(t). This result seems to be inconsistent with the taxonomic data, because our 
lifetime is measured by a log-scaled time. The apparent paradox is resolved by taking a 
small perturbation into account. In fact, we can show that a certain class of perturbations 
alters the expression of the scaled time r from r = log(t) to r = et while keeping the r 
dependence of statistical expressions unchanged [[□]]. e here measures the slowness of the 
time scale which is related to the magnitude of the perturbation. A similar change of a time 



scale is observed at a global bifurcation called a 'saddle connection' |L6 



In order to present more concrete discussion, we adapt [year] as the time unit of the 
population dynamics and set e = 10~ 6 as a typical value. Then, our result implies that (1) 
for species shorter lived than about 10 7 [year] extinction occurs randomly with respect to 
the age of the species at a constant rate of about 3 x 10 -7 [year -1 ], and (2) for species much 
longer lived than about 10 7 [year], older species have a smaller probability of extinction. The 
first result agrees with the taxonomic data of Van Valen ||. Also, there are some taxonomic 
data consistent with our second result PJTTf, though we have not checked the power law 
distribution for the lifetimes of longer lived species. According to our interpretation, the 
time scale of extinction is determined by slowness of small perturbations for population 
dynamics. We must consider what factors introduce the most relevant perturbations so as 
to estimate the value of e. Mutations of genes may be one candidate. It would be interesting 
to study the population dynamics of a system taking into account the effect of mutations 
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FIGURES 

FIG. 1. Time development of populations of N = 256 species. All populations within [0.01, 
0.99] are plotted every time unit. An initial condition was chosen randomly. An Euler method 
with a time mesh of 0.1 was adapted in the numerical simulation. 

FIG. 2. Plot of bursts in a species-r plane after 10000 bursts starting from an initial condition 
chosen randomly. 

FIG. 3. Distribution of intervals between two successive bursts for N = 256 species. The 
distribution is made from data of 40000 bursts after discarding the initial transient behavior. The 
distribution curve is well-approximated as P(t) = aexp(-ar), where a = 25. a would increase 
linearly with N. 

FIG. 4. Distribution of persistence time of positive hi. The distribution is made from data of 
80000 bursts after discarding the transient. The distribution curve has a power law tail for r > r c , 
while it is consistent with a Poisson distribution for r < r c . The cross-over time r c is about 10. 
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